Efficient and accurate calculation of exact exchange and RPA correlation energies in 
the Adiabatic-Connection Fluctuation-Dissipation theory 
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Recently there has been a renewed interest in the calculation of exact-exchange and RPA correla- 
tion energies for realistic systems. These quantities are main ingredients of the so-called EXX/RPA-I- 
scheme which has been shown to be a promising alternative approach to the standard LDA/GGA 
DFT for weakly bound systems where LDA and GGA perform poorly. In this paper, we present 
an efficient approach to compute the RPA correlation energy in the framework of the Adiabatic- 
Connection Fluctuation-Dissipation formalism. The method is based on the calculation of a rel- 
atively small number of eigenmodes of RPA dielectric matrix, efficiently computed by iterative 
density response calculations in the framework of Density Functional Perturbation Theory. We will 
also discuss a careful treatment of the integrable divergence in the exact-exchange energy calculation 
which alleviates the problem of its slow convergence with respect to Brillouin zone sampling. As 
an illustration of the method, we show the results of applications to bulk Si, Be dimer and atomic 
systems. 
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I. INTRODUCTION 

Developments during the last several decades have 
brought Density Functional Theory (DFT) to be the 
mcthod-of-choicc in many calculations of physical and 
chemical properties from first-principles. Despite its 
many spectacular successes, Kohn-Sham (KS) density- 
functional theory within the LDA and GGAs approxi- 
mations for the exchange-correlation (xc) functional still 
has some well-known drawbacks. This includes the fail- 
ure of standard DFT functionals in the description of 
van der Waals systems where long-range correlation ef- 
fects arc predominant. This failure manifests clearly in 
erratic or completely failed results when LDA/GGA DFT 
is applied to sparse matter, including layered structures, 
polymer and molecular crystals, or weakly-bound chem- 
ical systems. 

Of all the attempts to overcome these drawbacks, the 
approach based on the formally exact expression of xc en- 
ergy in term of linear response functions derived from the 
Adiabatic Connection Fluctuation-Dissipation (ACFD) 
theoreroi provides a promising way to develop a system- 
atic improvement. 

The method is in principle entirely parameter free- 
although an appropriate definition of the xc kernel is 
still needed-and a few difficult cases where standard 
DFT fails qualitatively have been described with satis- 
factory resultsi^i^'^ Practical applications in a few case- 
studies have shown that not only the description of van 
der Waals systems is improved but that the already sat- 
isfactory description of many molecules and solids are 
also not compromised.'^'^ The method is very compu- 
tationally demanding and this has prevented its direct 
application to complex systems. Most often it is lim- 
ited to a post-self-consistent correction where the accu- 



rate exchange-correlation energy is computed from the 
charge density obtained from a self-consistent calcula- 
tion performed with a more traditional xc functional. 
The ACFD formalism also serves as the starting point 
for further simplifications that have led to the develop- 
ment of an approximate vdW-DF functional by Langreth 
and coworkers which allows for the treatment of large 
systems and has recently been made fully self-consistent ii 

We will address here the computational difficulties 
involved in the full evaluation of ACFD exchange- 
correlation energies of realistic systems, presenting an 
efficient general formulation of the problem and its prac- 
tical implementation in the plane-wave pseudopotential 
(PW PP) formalism. 

Previous implementations of ACFD formulas rely on 
the construction of the full Kohn-Sham response function 
from the spectrum of the KS Hamiltonian. Basic defini- 
tions and some details are given in Sec. [TTl A more effi- 
cient scheme is then presented in Sec. Illll that avoids the 
full diagonalization of the KS Hamiltonian and the cum- 
bersome summation over occupied and unoccupied states 
needed in the traditional calculation of the response func- 
tions, as well as the computationally demanding opera- 
tions involved in the solution of the Dyson equation via 
exact algebra. These steps are replaced by the iterative 
calculation of a relatively small number of eigenmodes 
of the RPA dielectric matrix as explained in its general 
formulation in the first part of this section. Some tech- 
nical details of the implementation in the PW PP for- 
malism and a careful treatment of the integrable diver- 
gence which alleviates the problem of slow convergence 
with respect to Brillouin zone sampling in the calcula- 
tion of the exact-exchange energy are also discussed here. 
To demonstrate the efficiency of the implementation, we 
will present in Sec. II VI results of applications to selected 
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paradigmatic systems, namely bulk Si, weakly bound Be 
dimer, and a number of spherical atomic systems. 

II. THEORY 

We recall here some basic definitions in ACFD theory 
as well as some computational details of currently avail- 
able PW PP implementations. 

A. Exchange and correlation energy in ACFD 
theory 



an exchange-correlation kernel, ff^{iu), which is still un- 
known: 

=Xo+Xo[Awc + /r]XA- (5) 

In the above formula, frequency dependence is implied 
in each term (with the exception of the Coulomb kernel, 
Vc) and spatial coordinate dependence is implicit in the 
matrix notation. Despite the exactness of the expressions 
presented above, one needs to use an approximation for 
the xc-kernel in practical applications. When the xc- 
kernel is specified, the system of Eqs. ([3|) — ([5]) is closed 
and allows one to evaluate the correlation energy. 



In the ACFD formalism, exchange-correlation energy, 
Exci of an electronic system with density n(r) can be 
written in the formally exact form 



K-r = - ^ I dX I drdr'- ^ 
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r - r' 



X |y duxA(r,r';TO) + - r')n(r) j , (1) 

where x;^(r,r';?u) is the density response function at 
imaginary frequency, iu, of the system where Coulomb in- 
teraction between electrons is scaled by A, i.e. Ae^/|r— r'|, 
and the external potential is modified so that the elec- 
tronic density is the same as in the ground state of the 
physical system (A = 1). The exchange-correlation en- 
ergy can be furthermore separated into the KS exchange 
energy. Ex, and the correlation energy, E^. The former is 
the counterpart of Hartree-Fock exchange energy in the 
context of density-functional theory. It is given by the 
well-known expression in term of occupied KS orbitals (a 
non magnetic system is considered for simplicity and a 
factor of two accounts for spin degeneracy) 



drdr' 



ErV*(r)0.(r')| 



(2) 



The latter can be expressed in a compact form in terms 
of linear density responses, in matrix notation, as 



Er 
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duTr {vclxxiiu) - Xo{iu)]) . (3) 



In Eq. ([3]), Wc is the electron-electron interaction kernel, 
e^/|r — r'l, and Xo(*m) is the density response function of 
the non-interacting electron system and can be explicitly 
expressed in term of, occupied and empty, KS orbitals, 
4>i{v), KS eigenvalues, e,;, and occupation numbers, /,;, 

, ^^ 0*(r)0,(r)0*(rO^,;(rO 



B. Total energy in EXX/RPA+ scheme 

Experience has shown that none of the available ap- 
proximate xc-kernels gives a systematic improvement; 
different kernels perform better for different systems. Al- 
though the RPA-kernel, i.e. simply setting ff^"^- in Eq. ^ 
to zero, is exact for long-range correlation, it is a poor 
approximation for short-range correlation. When ap- 
plied to atomic or molecular systems, RPA correlation 
energy alone is very inaccurate. In fact, it was abandoned 
decades ago as a DFT approximate xc functional until a 
quite recent work by Kurth and Perdew^ showed that it 
is possible to correct the deficiency of RPA in a simple 
way by combining the full RPA with a local- or semilocal- 
density correction for short-range correlation. This leads 
to the so-called RPA 4- correlation energy in which the 
local-density correction for the short-range correlation is 
chosen in such a way that E^^^^ becomes exact in the 
limit of homogeneous electron gas (HEG) 



E. 



RPA+ 
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RPA 



LDA-RPA 
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LDA 



(6) 



In the above definition, £'LDA-rpa ^-j^^ local-density ap- 
proximation of RPA correlation energy which is exactly 
canceled out by Ef^^ in the limit of an HEG. 

In order to have exchange and correlation energies 
treated on the same footing, exchange energy must also 
be computed in the same formalism, i.e. from Eq. 

In practice this EXX/RPA-I- scheme is applied as a 
post-scf correction and the energy is calculated by first 
performing a standard LDA/GGA DFT calculation and 
then replacing the xc-energy at LDA level by exact- 
exchange and RPA+ correlation energies calculated from 
the LDA/GGA charge density. In spite of this, the 
EXX/RPA-h scheme has been shown to be a promising 
approach for the description of weakly bound systems 
where standard LDA/GGA performs poorlyi^i^ii^ 



C. Existing implementations in PW PP method 



where again a factor of two is present accounting for spin 
degeneracy. For A > 0, the interacting response func- 
tion, x\j is related to xo by a Dyson-like equation with 



Existing implementations^ of the ACFD formulas in 
plane-wave pseudopotential computer codes start by first 
diagonalizing the KS Hamiltonian for all the occupied 
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and unoccupied KS orbitals^ or at least a good part of 
the unoccupied states j^i^ so that the KS response func- 
tion xoiiu) can be calculated explicitly according to its 
definition in Eq. The calculation then proceeds by 
solving the Dyson- like equation for xa(*w) for a num- 
ber of values of the coupling constant, A, and of the 
imaginary frequency, iu, and then by integrating over 
these variables to ultimately obtain the correlation en- 
ergy. An obvious disadvantage of these implementations 
is that many unoccupied states need to be considered 
in order to get well-converged results. This forces one 
to solve the KS problem using full matrix diagonaliza- 
tion algorithms with unfavorable scaling instead of very 
efficient iterative-diagonalization techniques, commonly 
used to calculate the occupied KS orbitals in typical self- 
consistent calculations. Moreover the summation over 
valence and conduction bands and over all the k-points 
in the Brillouin zone for setting up the response function 
(Eq. (jH)) has been shown to be a very cumbersome op- 
eration which prevents the application of the method to 
large systems. 

III. EFFICIENT CALCULATION OF ACFD 
FORMULAS 

We describe now our method, first focusing on its gen- 
eral aspects, valid for any symmetry and basis set, and 
then specializing the discussion to a PW PP approach. 

A. General aspects and PW PP implementation 

Let us starts by defining the following generalized 
eigenvalue problems 

Xo\m) = a^v~^\w^), (7) 
Xx\z^) = bfv-^\z^), (8) 

where xo , Xa and Vc are linear response and Coulomb op- 
erators in matrix notation, {\wi),ai} and {\z^),b^} are 
eigenpairs (all these quantities depend implicitly on the 
imaginary frequency iu). These eigenvalues problems are 
well-defined since all the operators are symmetric and 
positive (vc) or negative (xo,Xa) definite. Once solu- 
tions of the generalized eigenvalue problems are avail- 
able, the traces of WcXa and VcXo in Eq. ([3]) are simply 
evaluated by summing up all their (relevant) eigenvalues. 
Working with these generalized eigenvalue problems has 
several advantages when RPA approximation is consid- 
ered: (i) the set of formally A-dependent eigenpotentials, 
{\z^)}, is actually A-independent and coincides with its 
non-interacting counterpart, {\wi)}, (ii) the correspond- 
ing eigenvalues are trivially related as 



which (iii) allows one to perform the A coupling constant 
integration in Eq. ^ analytically leading to the final 



expression 

Ec = 7r duY{a,{iu) + lii{l-a,{iu))}. (10) 

The spectrum of the eigenvalue problem in Eq. ([7]) 
is bounded from above by zero, which is also an eigen- 
value, corresponding to a constant eigenpotential. Note 
that the eigenvalues {ai} are closely related to those of 
the RPA dielectric matrix since eRPA = 1 — f^cXo, and 
the calculation of {ai} is very similar to the calculation 
of dielectric band structures introduced a few decades 
agOf^^iii In principle, the full spectra must be known in 
order to calculate the correlation energy from Eq. pU]) 
above. In practice, however, previous calculationsiSiiiii^ 
have shown that only a small number of eigenvalues of 
the RPA dielectric matrix significantly differs from unit. 
This translates in the fact that only a small fraction of 
the eigenvalues {a;} will be significantly different from 
zero and needs to be explicitly included in the correla- 
tion energy sum; all the rest being so close to zero that 
their contributions to Ec can be treated in a suitable ap- 
proximation or even simply neglected. 

Evaluation of low-lying eigenvalues of the non- 
interacting problem in Eq. ([7]) is done efhciently by it- 
erative diagonalization procedure starting from a set of 
trial eigenpotentials. The basic operation involved in the 
iterative diagonalization is the calculation of the non- 
interacting response to a trial potential, An = XoAi;. 
This is done resorting to the linear response techniques 
of Density Functional Perturbation Theory^H simply gen- 
eralized to imaginary frequency. 

Technically, for any given non-interacting perturtbing 
potential. Aw, the induced charge density variation is 

occ 

An(r) = 2Re^0*(r)A0,(r), (11) 

i 

where the sum runs over the occupied (valence) states 
and A0i(r) is the (conduction-band projected) variation 
of the single-particle state, (j)i{r), that can be obtained 
as the solution of the linear system: 

[Hks + aPv - {e^ + lu)] |A0,) = -{l-P,)Av\(j>,) , (12) 

where Py = I't'i) {4'i\ is the projector on the occu- 

pied manifold and a is a positive constant, larger than 
the valence bandwidth, ensuring that the linear system 
is not singular even in the limit for iu — *■ 0. The non- 
hcrmitian linear systems are solved by employing a fast 
and smoothly converging variant of the bi-conjugate gra- 
dient algorithmfi^ 

Overall the procedure amounts to solving the usual 
DFPT linear systemsi^ in which the ground state valence 
eigenvalues are shifted by the complex constant iu, as 
already explained by Mahartii long time ago for the case 
of atomic polarizabilities. Note however that, at variance 
with standard DFPT, no self-consistent cycle is needed to 
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obtain An as the response of the non-interacting system 
is being considered here. 

As in static DFPT calculations, only occupied Kohn- 
Sham orbitals and their linear response need to be cal- 
culated, instead of the full KS spectrum needed in other 
implementatiouj^i^ii^ 

The scheme describe above have been implemented 
in the Quantum ESPRESSO distribution^^ Very simi- 
lar ideas have also been recently reported by Wilson and 
coworkersii for the calculation of dominant eigenpoten- 
tials of static dielectric matrices. 

The description presented so far is to a large extent 
basis set independent. For periodic systems, matrix rep- 
resentations of the Kohn-Sham linear response functions 
satisfy the Bloch theorem and can be classified by a vec- 
tor q in the Brillouin zone. The generalized eigenvalue 
problem in Eq. ^ can therefore be solved separately for 
each q vector. As a consequence, the expression for the 
correlation energy per unit cell becomes 



Ti f 1 

Ec^— du — '^'^{a^{q,iu)+ln{^-a^{q,iu))}, 







(13) 

where the notation X]g=i indicates the average over 
the Brillouin zone and can be performed using the special 
point technique that takes advantage of the point group 
symmetry of the system. 

Special care is required in the limit of q since the 
leading matrix elements Xo^Iq ~* 0) and both go to 
zero as |qp making the eigenvalue problem ([7]) ill-defined. 
A simple solution to the problem is to use shifted grids 
of q vectors in the special-point summation. 



B. Exchange energy 

Although the main focus of this work is the calculation 
of correlation energy it is clear that an accurate evalu- 
ation of the exchange term is also necessary. A subtle 
problem, not widely recognized in the literature, related 
to the q ^ divergence of the Coulomb interaction is 
present in the calculation of exact-exchange energy in a 
plane-wave approach and special care must be payed to 
it. 



The exchange energy per unit cell for a periodic system 
is defined as 



-y. 



0*,(r)0k'.'(r)<^^,.(r')0k.(r') 



drdr' 



(14) 

where an insulating and non magnetic system is assumed 
for simplicity. Integrals and wave-function normaliza- 
tions are defined over the whole crystal volume, V = Nil 
(f2 being the unit cell volume), and the summations run 
over all occupied bands and all N k-points defined in the 
BZ by Born-von Karman boundary conditions. 

Introducing codensities, pk'.v'{r) = (r)0k,ti(r), 

with Fourier transform, p^'.v' (k — k' + G), the above ex- 
it, i- 

prcssion can be recast in reciprocal space as 



Er, 



(2-Sf/'^§ 



|q + Gp 



(15) 



where an auxiliary regular function j4(q -|~ G) = 
ir Ek J2v,v' IPk^^,'''" (q + G)P has been introduced. 

A direct integration of Eq. (fTS)) on the regular q-point 
grid, defined by the k — k' differences of the points used 
for the wave-function BZ sampling, is problematic due 
to integrable divergence that appears in the q + G ^ 
limit. The problem has been addressed several time in 
the literature and many related schemes have been pro- 
posed to treat it. Most eventually resort to a procedure, 
first proposed by Gygi and Baldereschiji^ where an eas- 
ily integrable term, that displays the same divergence, is 
subtracted from and separately added back to Eq. (|15p. 
Of the many proposed form a^^i^^i^° we adopt here one 
of the simplest that can be straightforwardly applied to 
any crystal structure: we simply subtract and add a term 
A(0)e-"li+'^lV|q + Gp from the integrand. The inte- 
gration of the resulting non divergent term can then be 
performed by a summation over a regular grid of points 
(either the full q-point grid defined from the k-point dif- 
ferences or, in order to reduce computational cost, on 
a subgrid of it), while the divergent integration can be 
performed explicitly. 

The final results is 



47re^ 



E 



'A{q + G) 



DA{0) 



(16) 



with 



D 



'g-a|q+G 



q,G 



G|2 



(2^ 



(17) 
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The symbol J2 the preceding formulas implies that 
the term q + G = is excluded. Notice that with our 
simple choice the expression for the divergence correc- 
tion, D, is just (the reciprocal-space part of) the Ewald 
sum for a single point-charge, periodically repeated ac- 
cording to the super-periodicity defined by the q-point 
grid used. As long as the parameter a is chosen so that 
the reciprocal space sum is converged the result is in- 
dependent of a and gives a correction that decay very 
slowly with the inverse of the cell linear dimension, L. 
This 1/L dependence of the singularity correction is for 
instance evident in Fig. 2 of Ref. [13, although the scal- 
ing law of the correction was not explicitly pointed out 
there. In general, the error incurred on by neglecting this 
correction would be so large that it is almost invariably 
included in the calculation. 

From a closer analysis of Eq. however, it is 

also clear that, even when the main effect of the di- 
vergence has been taken into account by the Gygi- 
Baldereschi procedure, another correction needs to be 
included to properly describe the spherically-averaged 



limit, ((limq^o 



^(q)-^(o) 



)), that cannot be calculated di- 



rectly for q = 0. 

It can be shown, see the Appendix, that this term is 
simply related to the gauge invariant spread of the occu- 
pied manifold, as can be defined in the theory of max- 
imally localized Wannier functions If this term is ne- 
glected (as to the best of our understanding it has always 
been neglected so far in the literature) an error propor- 
tional to this spread and inversely proportional to the 
cell volume and the number of q-points included in the 
BZ summation is made. This might be the reason for the 
reported slow convergence of the exchange energy with 
respect to BZ integration, even for simple systems such 
as bulk Silicon^ or Argoni^ 

We have therefore included in the calculation of the 
exchange energy an estimate of the limiting term, based 
on the assumption that the grid of q-points used for BZ 
integration is dense enough that a coarser grid, including 
only every second point in each direction would also be 
equally accurate. Since the limiting term contribute to 
the integral with different weight in the two grids one can 
estimate its value as: 



{(lim 



^(q) - ^(0) 



' ^(q + G) 
|q + G|2 



• ^(q + G) 
|q + G|2 



(18) 



r 



We have verified that the convergence with respect to 
BZ integration is generally improved. For instance for 
the exchange energy of solid Argon, calculated from LDA 
wave- functions, a regular grid of 6x6x6 points was suffi- 
cient for a convergence of the order of 0.1 niRy, while a 
much denser 12x12x12 grid was needed in Ref. [a- 



C. Computational cost 

To discuss the computational cost of our implemen- 
tation, let us denote by A'^, and A^c the number of va- 
lence and conduction bands, and by Np^^, and A'p^,^^ 
the number of plane waves used to represent wave func- 
tions and Kohn-Sham response functions, respectively. 
The basic operations in our implementation is the cal- 
culation of the linear density response via DFPT tech- 
nique. For a norm-conserving pseudopotential, the com- 
putational cost of a linear response calculation is essen- 
tially X {Ni, X Npw^ -t- Npt^^g). Since this calcula- 
tion is done repeatedly in the iterative diagonalization 
procedure, the total computational cost must be multi- 
plied by the number of eigenvalues (Neig) that we want 
to calculate and the number of iterations Nuer m the 
iterative diagonalization. While the latter is likely in- 
dependent of the system size, it is expected and in fact 
there is evidenceii that the former is proportional to it. 



Therefore the total scaling of our approach is propor- 
tional to Npuj^N^Neig which grows as the fourth power 
of the system size. In other implementations of ACFD 
formulas which are based on the evaluation of the full 
response matrix, the most time-consuming operation is 
the construction of the non-interacting Kohn-Sham re- 
sponse function whose computational cost is proportional 
to Np^^^N^Nc- This means that the computational cost 
of these approaches also scales with the fourth power of 
the system size. 

Note however that Npy^^, is typically smaller than 
Npwxo by an order of magnitude since xo relates to den- 
sity responses and perturbing potentials whose kinetic- 
energy cutoff is four times larger than the one needed 
for wave-functions in the case of norm conserving pseu- 
dopotentials and even more than that in the case of 
ultra-soft pseudopotentials. Also the number of cigcn- 
potcntials that needs to be computed is expected to be 
at least an order of magnitude smaller than the size of 
the response function (Np^^^) and the number of iter- 
ation Niter is unlikely to exceed 10, especially if in the 
imaginary frequency scan good starting trial eigenpoten- 
tials are taken from the previous frequency in the list. 
Our rough estimate shows therefore that the number of 
operations involved in our method is from 100 to 1000 
smaller than the one needed for other implementations. 
Although real CPU-time obviously depends on many de- 
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tails in the realization of each approach, we are confident 
that our approach allows us to implement ACFD formu- 
las in a very efficient way. 

Moreover, our approach also has significant advantages 
due to its iterative nature. Iterative solution of the gen- 
eralized eigenvalue problem ([8]) will converge very rapidly 
if the initial guess of the eigenpotentials are already close 
to the sought solutions. As mentioned earlier calculated 
solutions for a given imaginary frequency can serve as 
starting points for nearby frequencies. A similar behavior 
can be expected also for calculation of response functions 
in nearby q-vectors although it may be more convenient 
in this case to exploit the very easy parallelization of the 
q-vector summation involved in the correlation energy 
formula, Eq. (fT5|) . Finally, as already mentioned in the 
general formulation, the analysis of the response function 
in terms of eigenmodes allows us to analytically perform 
the coupling constant integration with a significant com- 
putational saving. 



IV. APPLICATION TO SELECTED SYSTEMS 
A. Bulk Si 

We have chosen bulk Si system as a testing ground for 
our implementation of ACFD theory since the compu- 
tational cost for the case of bulk Si is rather moderate, 
which is convenient for convergence checks, and there ex- 
ist several published data that we can compare with. As 
a first check, we have calculated the 20 topmost eigenval- 
ues of the dielectric matrix at q = (0, 0, 0.01 and com- 
pared our results with those reported in Ref. |l3 where the 
same quantities were calculated by the explicit diagonal- 
ization of the RPA dielectric matrix and by an iterative 
diagonalization procedure, named Projective Dielectric 
EigenPotential method (PDEP), similar to oursi^ The 
results reported in Table |T] show perfect agreement of the 
two methods among each other and with the explicit cal- 
culation. 

Let us now investigate the convergence of RPA correla- 
tion energy with respect to the number of eigenmodes in- 
cluded in the summation of Eq. (|13|) . To this purpose, we 
have used well-converged parameters for standard LDA 
calculation of silicon ground-state charge density to eval- 
uate the exact-exchange and RPA-correlation energies. 
The calculation was performed for the diamond structure 
in the fundamental face-centered cubic cell with lattice 
constant of 10.20 Bohr (corresponding to the theoretical 
equilibrium geometry at LDA level) using a regular grid 
of 64 k-point and a kinetic-energy cut-off of 20 Rydberg. 

We show in Fig. [1] the dependence of RPA correla- 
tion energy, Ej}^^, on the number of eigenvalues, N^ig, 
included in the summation. E^^^ is indeed a rapidly 
converging function of N^ig ; truncating the sum after in- 
clusion of 80 or 100 eigenvalues already ensures a con- 
vergence within a few tens of mRy. Also the summation 
over special q-points representing the integration in the 



Index 


RPAii 


PDEP^^ 


present work 


1 


14.7432 


14.7538 


14.7611 


2 


3.4231 


3.4237 


3.4238 


3 


3.3908 


3.3914 


3.3915 


4 


3.3908 


3.3914 


3.3915 


5 


3.3908 


3.3914 


3.3915 


6 


3.3908 


3.3914 


3.3915 


7 


3.3589 


3.3596 


3.3596 


8 


2.4910 


2.4925 


2.4925 


9 


2.4910 


2.4925 


2.4925 


10 


2.4910 


2.4925 


2.4925 


11 


2.4910 


2.4925 


2.4925 


12 


2.4905 


2.4920 


2.4920 


13 


2.4905 


2.4920 


2.4920 


14 


2.4716 


2.4721 


2.4721 


15 


2.1964 


2.1972 


2.1972 


16 


2.1960 


2.1968 


2.1968 


17 


2.1959 


2.1966 


2.1967 


18 


2.1959 


2.1966 


2.1967 


19 


2.1958 


2.1966 


2.1967 


20 


2.1958 


2.1966 


2.1967 



TABLE I: Top 20 eigenvaules of the dielectric matrix for an 
eight-atom silicon cubic cell calculated via the present itera- 
tive method and compared with he results reported in Ref. [l3 
for RPA and for projective dielectric eigenpotential method 

first Brillouin zone converges very rapidly; the correlation 
energy changes only by a few mRy when the number of 
special points increases from 2 to 6. 

Next we compare our accurately calculated exchange 
and RPA(-I-) correlation energies with the energies cal- 
culated within the pseudopotential approximation but 
using quantum Monte Carlo techniques as reported by 
Hood and co-workers in Ref. IH. Table [TTl shows clearly 




20 40 60 80 100 120 140 160 180 200 
Number of eigenvalues N . 

to eig 



FIG. 1: RPA correlation energy as function of the number of 
eigenvalues included in the summation in Eq. (|13|l for fee bulk 
silicon. The two curves are for different number of special q 
points in the Brillouin zone integration of the same equation. 
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ao (bohr) (eV) 



(eV) 



VMC 
DMC 

EXX/RPA 
EXX/RPA+ 



10.20 
10.26 

10.20 
10.26 



-29.15 
-29.15 

-29.11 
-28.98 

-29.11 
-28.98 



-3.58 ±0.01 
-4.08 ±0.08 

-6.12 
-6.11 

-4.24 
-4.23 



TABLE II: Exchange and correlation energies for bulk Si cal- 
culated by our method compared to the values calculated 
using QMC techniques by Hood and co-workers reported in 
Ref.il. 



that RPA alone gives indeed too negative values for the 
correlation energy. When a simple local-density correc- 
tion is added (RPA+), the presently calculated corre- 
lation energy becomes much closer to the Monte Carlo 
values. Obviously, the comparison would be more con- 
vincing had we used the same pseudopotential used by 
Hood and co-workers. Nonetheless, we believe that the 
good agreement between our calculations and the results 
obtained by a very different method shows that our re- 
sults are at least in the right regime. 

We have also calculated the total energy of bulk sili- 
con as a function of its lattice constant in EXX/RPA± 
scheme in order to determine the corresponding equi- 



0.05 




9.60 9.80 10.00 10.20 10.40 10.60 10.80 11.00 

a (a.u.) 

FIG. 2: The total energy (per unit cell) differences at different 
lattice constants and that at the experimental one of bulk 
silicon calculated using different number of eigenvalues values 
included in the summation in Eq. (|10|) in our implementation 
of EXX/RPA± scheme. The inset is the true total energy, 
i.e. the values at the experimental lattice constant is not 
subtracted. 



librium properties. Fig. [2] shows the total energy per 
unit cell for bulk silicon as a function of lattice con- 
stant calculated including different numbers N^ig {N^ig = 
50, 100, 150 and 200) of response fimction eigenvalues in 
the summation in Eq. pO)) . For each curve the value 
calculated at the experimental lattice constant is sub- 
tracted. The resulting (nearly) coincidence of the differ- 
ent curves confirms the expectation that energy differ- 
ences are rather insensitive to the number of eigenvalues 
included in the RPA sum. This can also be appreci- 
ated in the inset, where the unshifted EXX/RPA± to- 
tal energies are shown: as more eigenvalues are included 
in the sum energy vs volume curves rigidly shift down 
and although complete convergence is reached only when 
including about 200 eigenvalues already 100 or even 50 
eigenvalues basically return the same structural proper- 
ties. 

Table IIIII collect the predicted equilibrium lattice pa- 
rameter oo, bulk modulus B and pressure derivative of 
the bulk modulus B' as a function of the number of eigen- 
values Nf-ig used to evaluate RPA correlation energies. 
The changes in these quantities when increasing N^ig 
from 50 to 200 is very small, of the order of 0.1% for 
flo and B and 1% for B' . We can therefore conclude that 
only a relatively small number of eigenvalues are needed 
in the evaluation of ACFD RPA correlation energy to get 
a good description of many equilibrium properties of the 
system. 

Examining the predicted equilibrium lattice constants 
in Table [1111 it seems that a more accurate (and sophis- 
ticated) treatment of xc-energy in EXX/RPA± scheme 
slightly worsens the agreement with experimental data as 
compared with the LDA results. Nevertheless, as already 
pointed out in Ref. |3, there are several points that can 
affect the final results: (i) EXX/RPA+ scheme is applied 
in a non self-consistent way using LDA Kohn-Sham or- 
bitals, (ii) the pseudopotential itself has been generated 
within LDA, and (iii) RPA± is among the simplest pos- 
sible approximations to the xc-kernel within ACFD for- 
malism. In view of these shortcoming we believe that the 







ao (a.u.) 


B (GPa) 


B' 


50 




10.155 


99.5 


4.22 


100 




10.158 


99.4 


4.21 


150 




10.162 


99.3 


4.19 


200 




10.166 


99.1 


4.17 




LDA 


10.235 


92.5 


4.16 




Expt" 


10.26 


99.2 


4.15 



'^Experimental values are taken from Rcf . [25l 

TABLE III: Predicted equilibrium lattice parameter, ao, bulk 
modulus, B, and pressure derivative of the bulk modulus, 
B' , as function of the number of eigenvalues Ndg used to 
evaluate RPA± correlation energies. The corresponding LDA 
and experimental values are also shown for comparison. 
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present results for Silicon can be considered satisfactory. 

In spite of the improvement of numerical efficiency in 
our implementation, we also observed a slow convergence 
of RPA correlation energy with respect to the kinetic- 
energy cutoff as reported in other implementations.— 
Efficient extrapolation schemes that allow for the evalua- 
tion of correlation energies in the limit of infinite energy 
cutoff, already proposed for those implementations, can 
be easily adapted to ours. As for the convergence with 
respect to BZ sampling, we notice, as explicitly shown in 
Ref . [26I for the case of bulk Si, that the difference between 
RPA correlation energies calculated using different num- 
ber of k-points (for charge density) and q-points (for Ec) 
is a well-behaved function of the kinetic-energy cutoff, 
beyond a certain (not very large) value. This suggests 
an extrapolation scheme, similar to the one proposed in 
Ref. 13, that might become useful for the calculation of 
more complex systems for which even the present very ef- 
ficient implementation would be too demanding to reach 
complete convergence directly. The extrapolation proce- 
dure would be as follows: First a coarse grid of k- and 
q-point could be used for the calculation of RPA correla- 
tion energies at different kinetic-energy cutoffs. Second, 
the corresponding correlation energy in the infinite cutoff 
limit could be obtained by extrapolating the results ob- 
tained at finite cutoffs. Finally, the errors due to coarse 
k- and q-point sampling of the BZ could be corrected by 
using finer grids evaluated at small kinetic-energy cutoff 
(whose safe value could be estimated from the conver- 
gence behavior of the correlation energy computed with 
the coarser grids.) 



B. Be2 dimer 

Beryllium dimer is a paradigmatic example of the fail- 
ure of LDA/GGA DFT in the description of weakly 
bound systems. Previous studiea^ii^ have shown con- 
siderable discrepancies between LDA/GGA and experi- 
mental results for binding energy, bond length and vibra- 
tional frequencies. While the errors of LDA/GGA bond 
lengths are in fact quite small (less than 2%), the vibra- 
tional frequency is largely overestimated. The most se- 
vere discrepancy refers to the binding energy: both LDA 
and GGA approximations overestimate the experimen- 
tal value by at least a factor of 4. Fuchs and Gonzo^ 
have recently shown that EXX/RPA-I- improves signifi- 
cantly the description of Be2 dimer over standard DFT. 
However, the error-bar of the RPA correlation energy re- 
ported in that study is still as large as 30 meV which 
is of the same order of the experimental binding energy 
and the reported theoretical value. Our efficient imple- 
mentation of ACFD formulas allows us to reach a better 
congergence, and thus to provide a better assessment of 
the performance of EXX/RPA-I- formalism for this sys- 
tem. To this end, we have carefully checked the con- 
vergence with respect to all parameters involved in our 
calculations, namely the kinetic-energy cutoff, the size 
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FIG. 3: Exact-exchange and RPA correlation energies of Be2 
as a function of Be-Be distance (dse-Se) as calculated in a 
simple cubic cell of side length of 22 (square) and 30 (circle) 
Bohr. Dotted lines are simply drawn as a guide to the eye. 



of the supercell used to simulate an isolated system us- 
ing a periodic plane-wave approach, and the number of 
eigenvalues included in the evaluation of RPA correlation 
energy. 

The same norm-conserving pseudopotential used in 
Ref. 0, generated within exact-exchange Kohn-Sham for- 
malism, has been used here. It is a rather soft pseu- 
dopotential and increasing the plane-wave kinetic-energy 
cutoff from 20 to 30 Ry makes the LDA total energy, the 
exact-exchange and RPA correlation energies change less 
than 0.5 mRy. This good convergence is consistent with 
the use in Ref. of a plane- wave kinetic-energy cutoff of 
25 Ry. 

Convergence with respect to the size of supercell is 
more delicate. Figure [3] shows exact-exchange and RPA 
correlation energies as a function of Be-Be distances cal- 
culated using supercells of 22 and 30 Bohr. While RPA 
correlation energy varies by only a fraction of a mRy 
(top panel) , the lower panel shows a significant variation 
of several mRy for the exact-exchange energy. As dis- 
cussed in Sec. IIII Bl slow convergence of exact-exchange 
energy with respect to supercell size-or the density of 
the BZ sampling for an extended system-can result if, 
once the integrable divergence is eliminated by the Gygi- 
Baldcrcschi proccdureji^ the residual q = term is not 
estimated correctly. This is shown in Figure |4] where 
the exact exchange energy of Be2 molecule is shown as a 
function of the inverse supercell volume. A large error, 
proportional to the inverse supercell volume, is present 
when the residual q = term is simply neglected, while 
a much better convergence with system size is obtained 
when it is estimated according to the recipe described 
in Sec. IIIIBI While with a simple cubic supercell of 22 
Bohr side the calculated exchange energy still has an er- 
ror of a couple of mRy, we are fully confident that using a 
supercell of 30 Bohr ensures a convergence of the exact- 
exchange energy within a fraction of mRy. One might 
expect some degree of error cancellation in energy dif- 
ferences when the residual q = terms are neglected in 
the exact-exchange calculations of Be dimer and of the 
separated Be atoms. Our explicit verification has shown 
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3 

cd 



with correction 
without correction 



^ 2x10 4x10 6x10 8x10 1x10 

1 /volume 

FIG. 4: Exchange energy of Be2 as a function of the supercell 
volume. Both the calculation neglecting (filled circles) and 
including (filled squares) an estimate of the q = term in the 
sum are reported. Numbers close to the symbols correspond 
to the supercell linear dimension in Bohr. 



that a variation of the order of mRy is stiU present in 
this difference when the supercell size is increased from 
22 to 30 bohr. 

Finally, convergence of RPA correlation energy with re- 
spect to the number of eigenvalues included in the ACFD 
summation in Eq. (fTU|) is shown in Fig. [5] where RPA cor- 
relation energies of Be2 calculated using 120, 180, and 220 
eigenvalues are plotted as a function of Be-Be distances. 
By including up to 220 eigenvalues, a convergence within 
0.5 mRy is obtained for the absolute value of correla- 
tion energy. On the basis of the nearly parallel behavior 
of the curves reported in Fig. [5] for different number of 
included eigenvalues we can anticipate a much faster con- 
vergence for energy differences, which actually determine 
the equilibrium properties of the dimer. 

Table HV] collects predicted binding energy, bond 
length, and vibrational frequency of Be2 calculated in- 
cluding different numbers of eigenvalues in the ACFD 
evaluation of RPA correlation energy. As expected, all 
these quantities are only very slightly changed when the 



E5(eV) 



do (bohr) 



.(cm-i) 





RPA 


RPA-I- 


RPA 


RPA + 


RPA 


RPA + 


60 


-0.0667 


-0.0377 


4.516 


4.553 


296.1 


298.5 


120 


-0.0657 


-0.0368 


4.521 


4.558 


296.6 


298.7 


180 


-0.0655 


-0.0367 


4.523 


4.560 


296.3 


299.2 


220 


-0.0654 


-0.0365 


4.524 


4.561 


297.1 


298.5 




FIG. 5; RPA correlation energies of Be2 at difi'erent Be-Be 
distances (dse-Se). The curves are for different number of 
eigenvalues Neig included for evaluation of Eq. (|10|l : A'' = 120 
(circle), N = 180 (square), and N = 220 (diamond). The lines 
are simply drawn as a guide. Inset: RPA correlation energies 
of Be2 at the Be-Be distance (dse-Se) of 4.56 Bohr placed 
in a simple cubic supercell with the size length of 22 (dashed 
line) and 30 Bohr (solid line). 



number of eigenvalues is increased from 60 to a value 
as large as 220. Even the molecular binding energy ap- 
pears to have a convergence of the order of the meV when 
only 60 eigenvalues are included in the ACDF sum. The 
result that properties related to energy differences are 
rather insensitive to the number of cigcnpotcntials used 
in evaluating RPA correlation energy, not only in the 
"standard" case of silicon (section [IV A[) . but also in the 
much more delicate case of Be2 dimer is very promis- 
ing for the future application of the present method to 
complex realistic systems. 

When compared with the reference values reported in 
Ref.i, our binding energies, both in RPA and RPA+, arc 
systematically smaller, although still in agreement within 
the large error bar of the previous calculation. Moreover, 
had we calculated the binding energy of Be2 without the 
q = correction described above, as done in the earlier 
calculation^ a much better matching with the ACFD 
calculation in Ref. would have been obtained. 

Coming to the comparison with the experimental re- 
sults it is clear that while EXX/RPA-f scheme definitely 
improves the poor performance of LDA and GGA for 
weakly bound systems like Be2, the results may not be 
as good as suggested in Ref. The performance of 
RPA-f- and other ACFD-based schemes to describe re- 
alistic weakly bound systems needs to be more system- 
atically investigated. 



Ref. -0.08(3) -0.06(3) 4.55 4.59 311 298 
Expt. -0.098 4.63 275.8 



TABLE IV: Equihbrium properties of Be2 in EXX/RPA(+) 
scheme with different numbers of eigenvalues included in the 
calculation of RPA correlation energy. 



C. Spherical atomic systems 

Benchmark results of RPA correlation energies for a 
number of spherical atoms calculated by constructing the 
full response function from the spectrum of Kohn-Sham 
Hamiltonian have been recently reported in Ref. [3^ . 
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TABLE V: Full RPA and RPA+ correlation energy (in Rydberg atomic units) of spherical atoms compared to the reference 
and exact values. The reference data were calculated from EXX-only (i.e. exact-exchange and no correlation) KS orbitals in a 
different implementation.— 



Atom 


^cxpt 


















^EXX 


Ref. 32 




pEXX 


Ref. 32 


He 


-0.084" 


-0.229 


-0.168 


-0.167 


-0.166 


-0.096 


-0.094 


0.094 


Be 


-0.190" 


-0.447 


-0.373 


-0.367 


-0.358 


-0.230 


-0.224 


0.216 


Ne 


-0.786" 


-1.474 


-1.216 


-1.195 


-1.194 


-0.821 


-0.800 


0.800 


Ar 


-1.463" 


-2.842 


-2.221 


-2.206 


-2.202 


-1.503 


-1.487 


1.482 


Kr 


-4.15' 


-6.533 


-5.226 


-5.192 


n/a 


-3.736 


-3.702 


n/a 


Xe 


-6.86' 


-10.358 


-8.312 


-8.278 


n/a 


-6.049 


-6.016 


n/a 



" Experimental values quoted in Ref. 

^ Differenee of "exaet" and HF total energies reported in Ref. 1311 . 



For spherically symmetric systems, ground state Kohn- 
Sham orbitals arc classified by their principal quantum 
number, n, and by angular momentum numbers /, m. 
The KS equations can be solved numerically on a ra- 
dial grid, within a given approximate, LDA or GGA, 
functional. Similarly, the non interacting and interact- 
ing response functions are block-diagonal with respect to 
angular momentum I, and {21 + l)-fold degenerate with 
respect to m. Thus, contributions to the ACFD formula 
from different angular momenta can be calculated inde- 
pendently and added up. For each angular momentum 
the calculation proceeds as follows, (i) A trial poten- 
tial is selected and the corresponding linear density re- 
sponse is calculated by solving the modified Sternheimer 
equationii^ A single iteration, with no self-consistency, 
is needed since the non-interacting response-function is 
studied, (ii) The calculated density response is orthogo- 
nalized. with overlap matrix Wc, with respect to any previ- 
ously computed one. The generating potential is accord- 
ingly transformed, exploiting the linearity of xo- By cal- 




Number of eigenvalues 



FIG. 6: The dependence of RPA correlation energy on the 
number of eigenvalues for Xenon atom. The curves show the 
cumulative contribution to RPA correlation energy for some 
of the lowest angular momentum numbers (higher I are not 
plotted). Including the first 15 eigenvalues is enough to en- 
sure a convergence within ImRy, which is also used as the 
threshold for convergence with respect to angular number I. 



culating the Hartree potential from the resulting density 
a new trial perturbation is generated, (iii) A matrix rep- 
resentation of Xo on the hence generated trial-potential 
basis is built, and diagonalized to get the eigenvalues a^'s. 
This three-step process is repeated until convergence in 
the sum over eigenvalues in Eq. ([TU]) is reached. The 
same calculation procedure is then repeated at different 
values of angular momentum, /, and imaginary frequency, 

Let us investigate the convergence of the RPA Ec with 
respect to the number of eigenvalues of the generalized 
eigenvalue problem of Eq. ([7]). Fig. [S] shows the depen- 
dence of RPA correlation energy-separated in the dif- 
ferent angular momentum contributions-on the number 
of eigenvalues included in the sum for Xenon atom, the 
heaviest atom considered in this work (56 electrons). Ba- 
sis set convergence is carefully checked and, for the case 
of Xe, a basis-set size of 25 trial potentials is enough 
for the desired accuracy. It is clearly seen that the cor- 
relation energy converges quite rapidly; including up to 
15 eigenvalues is enough in order to convergence the to- 
tal correlation energy within ImRy. These calculations 
therefore confirm explicitly, also for the case of spherical 
atoms, our expectation that RPA correlation energy can 
be obtained from only a small number of eigenmodes of 
the non-interacting response-function of the system. 

Table El shows the full RPA and RPA-f correlation en- 
ergies calculated with our method for a number of spher- 
ical atoms whose ground state densities have been gen- 
erated from EXX-only2£ and standard LDA functionals. 
Experimental^ or accurate theoretical^ values are also 
shown for comparison, together with reference RPA and 
RPA+ correlation energies calculated recently^^ within 
a different implementation of the ACFD approach, from 
EXX-only charge densities. All our calculated values in 
table |V] are converged within a few mRy. The results 
slightly depend on the quality of the ground-state den- 
sity, i.e. on the different approximate xc-functionals used 
in the self-consistent ground state calculation. Focusing 
on the results obtained starting from EXX-only charge 
densities, our calculated values for the full RPA correla- 
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tion energy (the fifth column) agree well with the refer- 
ence data (sixth column) within the error bar (with the 
exception of Be case); a similar agreement is found for 
RPA-I- correlation energies (eighth and ninth columns of 
table |V|. The small residual differences between the val- 
ues obtained in the present and in the reference calcula- 
tions may probably be attributed to some slight residual 
difference in the electronic densities used as input. 

Not surprisingly, RPA correlation energies alone 
largely overestimate the exact values. When combined 
with a local-density corrections to form the RPA-I- ap- 
proximation the correlation energy compare very favor- 
ably with the exact values thus giving support to the 
validity of RPA-I- scheme. 



V. CONCLUSION 

In this work, we have proposed an efficient method 
for the calculation of RPA correlation energies in 
the adiabatic-coupling fluctuation-dissipation formalism. 
Our approach involves the evaluation of a relatively small 
number of eigenvalues of the non-interacting response 
function, obtained combining concepts from density func- 
tional perturbation theory with iterative diagonahzation 
techniques. General strategies as well as technical de- 
tails of the method as implemented in the plane-wave 
pseudopotential Quantum-ESPRESSO distribution have 
been discussed to some extent. 

We have applied the method to study a few systems, 
representative of bulk solids, weakly bound molecules and 
atoms. While the study of bulk silicon crystal helps val- 
idate the implementation, our study of Beryllium dimcr, 
thanks to its improved numerical accuracy with respect 
to previous studies, allows us to gain a clearer picture 
of the performance of EXX/RPA+ scheme in describ- 
ing weakly bound systems. Our calculation confirms the 
important improvements of EXX/RPA+ with respect to 
LDA or GGA but also shows that its performance in 
delicate cases can be less impressive than previously con- 
cluded. The good agreement of RPA+ correlation energy 
with experimental or accurate theoretical results for a 
few spherical atoms lend anyhow support to the quality 
of RPA-f correlation energies. 

We are confident that the possibility of a careful con- 
trol of the numerical accuracy in AGED calculation, re- 
sulting from the improved numerical efficiency of our 
method, will turn out to be very useful in the needed 
analysis of the performance of such a sophisticated den- 



sity functional in the study of many other realistic sys- 
tems. 
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APPENDIX 

In this appendix we demonstrate that the spherically- 
averaged limit, ((limq^o ^^^^^r^^)), is simply related to 
the gauge invariant spread of the occupied KS manifold, 
il/ , defined in the theory of maximally localized Wannier 
functions^^ as 



E 



(07i|r2|0?i) - ^ |(0?i|r|RTO)|^ 



R.i 



where {|Rn)} is (any) set of Wannier functions that de- 
scribes the occupied manifold. 

Starting from the explicit expression for A(q) in terms 
of the occupied Bloch states 

^(q) = ^EE<'^k,.|e-^'i'-|0k-q..')(</'k-q,„'|e+'''"'|0k..), 

k v,v' 

let us expand them in terms of the Wannier functions, 
[./-k,.) =5]C/.,„(k)^e"'^|Rn), 

n R 

where ?7„.„(k) are general unitary transformations in the 
set of occupied bands at point k in the BZ. Inserting this 
definition in the expression for ^(q), this becomes, after 
some straightforward manipulation, 

A(q) = ^^(Rn|e-'i''|0"i)(Om|e+"i'"|R-")- 

R rt.ni 

Expanding now the exponential factors in powers of q 
and taking the spherically averaged limit for q — > one 
easily obtains the desired result: 

q^o q^ 3 
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